# R script for replicating results in appendix section "Further Results"
# Zhai & Garside 2022
# original device: MacPro 13, R 4.1.2
# recommended working directory: Desktop

setwd("~/Desktop") # default wkd
rm(list = ls())

# pkgs --------------------------------------------------------------------

if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse, magrittr, broom, sf, fixest, estimatr, ggthemes, ggpubr, texreg)

# helpers -----------------------------------------------------------------

source("R_00_helpers.R")

# setups ------------------------------------------------------------------

theme.default <- theme_get()
theme_set(theme_minimal(base_size = 14))

IDs.models <- c("Base", "Full")
IDs.measures <- c("Primary", "Secondary")

# figure 8 ----------------------------------------------------------------

df.main <- read.csv("data_vote_main.csv", row.names = 1) %>% 
  group_by(land_name) %>% 
  mutate(flooded_land = max(flooded, na.rm = TRUE)) %>% 
  ungroup() %>% 
  filter(flooded_land == 1) 

xx.h1 <- "+pop_per_sqkm+net_out_per_thousand+old_share+income_mean+unemployed_rate+land_set_pct+land_agri_pct+factor(incumbent)+distance*factor(date)"
fm.h1 <- "v_green_pct~factor(flooded)*factor(date)"
fm.h1.x <- paste0(fm.h1, xx.h1)
fm2.h1 <- gsub("flooded", "flooded2", fm.h1)
fm2.h1.x <- paste0(fm2.h1, xx.h1)

df.h1 <- map(
  list(fm.h1, fm.h1.x, fm2.h1, fm2.h1.x),
  ~fit_twfe(df = df.main, fm = .x, fe = "gen_id", cl = "gen_id")) %>% 
  map(~{coeftable(.x) %>% as.data.frame}) %>% 
  map(~{rownames_to_column(.x, var = "term") %>% select(term, mean = Estimate, se = `Std. Error`)}) %>% 
  map(~{filter(.x, term %in% c("factor(flooded)1:factor(date)2021-09-26", "factor(flooded2)1:factor(date)2021-09-26")) %>% 
      mutate(term = case_when(grepl("flooded2", term) ~ IDs.measures[2], TRUE~IDs.measures[1]))}) %>% 
  `names<-`(rep(IDs.models,2)) %>% 
  bind_rows(.id = "model") %>% 
  mutate(lwr = mean-1.96*se, upr = mean+1.96*se)

ggplot(df.h1) +
  geom_pointrange(aes(x="bhat", y=mean, ymin=lwr, ymax=upr, color=model, shape = model), position = position_dodge(0.7)) +
  facet_grid(rows = vars(term)) +
  geom_hline(yintercept = 0, color = "darkred", lty = "dotted") +
  coord_flip() +
  scale_color_stata() +
  labs(x = NULL, y = NULL, color = "Model", shape = "Model", title = "Direct Effect") +
  scale_y_continuous(limits = c(-0.005, 0.03), breaks = scales::pretty_breaks(n=5)) +
  theme(axis.ticks.y = element_blank(), axis.text.y = element_blank())
ggsave("~/Desktop/appendix_figure8.pdf", width = 7, height = 4)  

# table 4 -----------------------------------------------------------------

df.main.h1a <- 
  df.main %>% 
  mutate(nsevere = as.integer(flooded==1 & severe==0),
         nsevere2 = as.integer(flooded2==1 & severe2==0))

fm.h1a <- "v_green_pct~factor(severe)*factor(date)+factor(nsevere)*factor(date)"
fm.h1a.x <- paste0(fm.h1a, xx.h1)
fm2.h1a <- gsub("severe", "severe2", fm.h1a)
fm2.h1a.x <- paste0(fm2.h1a, xx.h1)

tab4 <- map(
  list(fm.h1a, fm.h1a.x, fm2.h1a, fm2.h1a.x),
  ~fit_twfe(df = df.main.h1a, fm = .x, fe = "gen_id", cl = "gen_id"))
screenreg(tab4, digits = 3)

# tables 5 & 6 ------------------------------------------------------------

## data ----
id.s2 <- 
  df.main %>% 
  filter(land_name %in% c("Rheinland-Pfalz","Nordrhein-Westfalen")) %>% 
  select(gen_id, land_name) %>% 
  distinct()
df.satellite <- st_read("data_flood_satellite/merged-floods.shp")
df.flood.s2 <- st_read("data_flood_govt/gemeinde-flood.shp") %>% 
  rename_with(tolower) %>% 
  filter(bez!="Gemeindefreies Gebiet") %>% 
  select(ags, gen, flooded, severe) %>% 
  mutate(ags = as.numeric(ags)) %>% 
  mutate(gen_id = case_when(
    gen == "Berlin" ~ 11,
    gen == "Hamburg" ~ 2,
    (ags %% 1e3 == 0) ~ ags/1e3,
    TRUE ~ ags
  )) %>% 
  transmute(gen_id = as.integer(gen_id), gen_name = gen,
            flooded = flooded, severe = severe) %>% 
  distinct() %>% 
  filter(gen_id %in% id.s2$gen_id)

## x = area ----
sf::sf_use_s2(FALSE)
area.intersect <- 
  st_intersection(st_make_valid(df.flood.s2$geometry), st_make_valid(df.satellite$geometry))  %>% 
  st_area()
area.all <- 
  st_area(st_make_valid(df.flood.s2$geometry)) %>% 
  setNames(df.flood.s2$gen_id)
id.intersect <- 
  st_intersects(st_make_valid(df.flood.s2$geometry), st_make_valid(df.satellite$geometry) ) %>% 
  `names<-`(df.flood.s2$gen_id) %>% 
  .[which(sapply(., length)>0)] 
df.intersect <- 
  unlist(id.intersect, use.names = FALSE) %>% 
  setNames(., rep(names(id.intersect), lengths(id.intersect))) %>% 
  cbind.data.frame(gen_id = names(.), sta_id = ., area = area.intersect) %>% 
  group_by(gen_id) %>% 
  summarise(area = sum(area)) %>% 
  ungroup() %>% 
  mutate(gen_id = as.integer(gen_id))
df.intersect.pct <- 
  area.all %>% 
  {.[names(id.intersect)]} %>% 
  data.frame() %>% 
  rename(area_all = 1) %>% 
  rownames_to_column(var = "gen_id") %>% 
  mutate(gen_id = as.integer(gen_id)) %>% 
  left_join(df.intersect, ., by = "gen_id") %>% 
  mutate(area_pct = area/area_all)
df.s2 <- 
  df.intersect.pct %>% 
  mutate(area_pct = as.numeric(area_pct)) %>% 
  mutate(area_pct_q4 = case_when( 
    area_pct < quantile(area_pct, 0.25) ~ 1L,
    area_pct >= quantile(area_pct, 0.25) & area_pct < quantile(area_pct, 0.5) ~ 2L,
    area_pct >= quantile(area_pct, 0.5) & area_pct < quantile(area_pct, 0.75) ~ 3L,
    TRUE ~ 4L
  )) %>% 
  mutate(area_pct1 = case_when( 
    area_pct < quantile(area_pct, 0.05) ~ quantile(area_pct, 0.05),
    area_pct > quantile(area_pct, 0.95) ~ quantile(area_pct, 0.95),
    TRUE ~ area_pct
  )) %>% 
  left_join(df.main) 

## models ----
tab5 <- fit_twfe(df = df.s2, 
                 fml = "v_green_pct~I(area_pct1^2)*factor(date)+area_pct1*factor(date)+pop_per_sqkm+net_out_per_thousand+old_share+income_mean+unemployed_rate+land_set_pct+land_agri_pct+factor(incumbent)+distance*factor(date)",
                 fe = "gen_id", cl = "gen_id")
screenreg(tab5, digits = 3)

tab6 <- fit_twfe(df = df.s2, 
                 fml = "v_green_pct~factor(area_pct_q4)*factor(date)+pop_per_sqkm+net_out_per_thousand+old_share+income_mean+unemployed_rate+land_set_pct+land_agri_pct+factor(incumbent)+distance*factor(date)",
                 fe = "gen_id", cl = "gen_id")
screenreg(tab6, digits = 3)

# table 7 & figure 9 ------------------------------------------------------

## data ----
df.main0 <- read.csv("data_vote_main.csv", row.names = 1) %>% 
  group_by(land_name) %>% 
  mutate(flooded_land = max(flooded, na.rm = TRUE)) %>% 
  ungroup()
df.spatmat <- readRDS("data_spatmat.RDS") 
df.spatmat.knn <- 
  df.spatmat$knn %>% 
  as.data.frame() %>% 
  mutate(gen_id = as.integer(colnames(.)))
var.h1b <- 
  df.main0 %>% 
  filter(date == max(date)) %>% 
  select(gen_id, land_name, flooded, flooded2) %>% 
  left_join(df.spatmat.knn, by = "gen_id")
d.h1b <- as.matrix(var.h1b[,names(var.h1b) %in% df.spatmat.knn$gen_id]) %*% var.h1b$flooded 
d2.h1b <- as.matrix(var.h1b[,names(var.h1b) %in% df.spatmat.knn$gen_id]) %*% var.h1b$flooded2
df.h1b <- 
  data.frame(
    gen_id = rep(var.h1b$gen_id, 2),
    flooded_w = rep(d.h1b, 2),
    flooded2_w = rep(d2.h1b, 2),
    date = c(rep("2017-09-24", length(d2.h1b)), rep("2021-09-26", length(d2.h1b)))
  ) %>% 
  left_join(df.main0)
df.h1b.s4 <- df.h1b %>% 
  filter(land_name %in% df.main$land_name)

## models ----
fm.h1b <- "v_green_pct~factor(flooded2)*factor(date)+flooded2_w*factor(date)"
fm.h1b.x <- paste0(fm.h1b, xx.h1)

## table 7 ----
tab7 <- map(
  list(fm.h1b, fm.h1b.x),
  ~fit_twfe(df = df.h1b.s4, fm = .x, fe = "gen_id", cl = "gen_id"))
screenreg(tab7, digits = 3)

## figure 9 ----
df.fig9 <- 
  tab7 %>% 
  map(~{coeftable(.x) %>% as.data.frame}) %>% 
  map(~{rownames_to_column(.x, var = "term") %>% select(term, mean = Estimate, se = `Std. Error`)}) %>% 
  map(~{filter(.x, term %in% c("factor(date)2021-09-26:factor(nsevere)1","factor(date)2021-09-26:factor(nsevere2)1",
                               "factor(severe)1:factor(date)2021-09-26",
                               "factor(severe2)1:factor(date)2021-09-26",
                               "factor(date)2021-09-26:factor(nsevere)1",
                               "factor(flooded2)1:factor(date)2021-09-26",
                               "factor(date)2021-09-26:flooded2_w")) %>% 
      mutate(term = case_when(
        grepl("flooded2_w", term) ~ "Neighbours",
        grepl("(flooded2|nsevere2)", term) ~ IDs.measures[2], 
        grepl("(flooded|nsevere)", term) ~ IDs.measures[1],
        grepl("severe2", term) ~ "Secondary+",
        grepl("severe", term) ~ "Primary+"
      ))}) %>% 
  bind_rows() %>% 
  add_column(model = rep(IDs.models, each=2)) %>% 
  mutate(lwr = mean-1.96*se, upr = mean+1.96*se) %>% 
  mutate(term4 = factor(term, levels = c("Secondary", "Neighbours"), labels = c("Direct (Secondary)", "Indirect")))
ggplot(df.fig9) +
  geom_pointrange(aes(x="bhat", y=mean, ymin=lwr, ymax=upr, color=model, shape = model), position = position_dodge(0.7)) +
  facet_grid(rows = vars(term4)) +
  geom_hline(yintercept = 0, color = "darkred", lty = "dotted") +
  coord_flip() +
  scale_color_stata() +
  labs(x = NULL, y = NULL, color = "Model", shape = "Model", title = "Indirect Effect") +
  scale_y_continuous(limits = c(-0.05, 0.05), breaks = scales::pretty_breaks(n=5)) +
  theme(axis.ticks.y = element_blank(), axis.text.y = element_blank())
ggsave("~/Desktop/appendix_figure9.pdf", width = 7, height = 4)  

# tables 8 & 9 ------------------------------------------------------------

## data ----
df.pcpt0 <- read.csv("data_gles_wknr.csv", row.names = 1)
df.pcpt <- 
  df.pcpt0 %>% 
  select(wk_id, year, climate_position_pre = climate_position, climate_salience_pre = climate_salience) %>% 
  filter(year==min(year)) %>% 
  select(-year) %>% 
  left_join(df.pcpt0, ., by = "wk_id") 
df.pcpt.s4 <- 
  df.pcpt %>% 
  group_by(land_id) %>% 
  mutate(flooded_land = max(flooded, na.rm = TRUE)) %>% 
  ungroup() %>% 
  filter(flooded_land == 1) 

## models ----
xx.pcpt <- "+pop_per_sqkm+net_out_per_thousand+old_share+income_mean+unemployed_rate+land_set_pct+land_agri_pct+factor(incumbent)+distance*factor(year)"
fm.pcpt.a <- "v_green_pct~factor(flooded)*factor(year)*climate_position_pre"
fm.pcpt.b <- "v_green_pct~factor(flooded)*factor(year)*climate_salience_pre"
fm.pcpt.x.a <- paste0(fm.pcpt.a, xx.pcpt)
fm.pcpt.x.b <- paste0(fm.pcpt.b, xx.pcpt)
fm2.pcpt.a <- gsub("flooded", "flooded2", fm.pcpt.a)
fm2.pcpt.b <- gsub("flooded", "flooded2", fm.pcpt.b)
fm2.pcpt.x.a <- paste0(fm2.pcpt.a, xx.pcpt)
fm2.pcpt.x.b <- paste0(fm2.pcpt.b, xx.pcpt)

## table 8 ----
tab8 <- map(
  list(fm.pcpt.a, fm.pcpt.x.a, fm2.pcpt.a, fm2.pcpt.x.a), 
  ~fit_twfe(df = df.pcpt.s4, fm = .x, fe = "wk_id", cl = "wk_id"))
screenreg(tab8, digits = 3)

## table 9 ----
tab9 <- map(
  list(fm.pcpt.b, fm.pcpt.x.b, fm2.pcpt.b, fm2.pcpt.x.b),
  ~fit_twfe(df = df.pcpt.s4, fm = .x, fe = "wk_id", cl = "wk_id"))
screenreg(tab9, digits = 3)

# figures 10 & 11 ---------------------------------------------------------

df.fig10 <- 
  df.main %>% 
  group_by(date, flooded, land_name) %>% 
  summarise(v_green_pct_m = lm_robust(v_green_pct~1)$coefficients,
            v_green_pct_lwr = lm_robust(v_green_pct~1)$conf.low,
            v_green_pct_upr = lm_robust(v_green_pct~1)$conf.high) %>% 
  ungroup() 
ggplot(df.fig10, aes(factor(date), v_green_pct_m, color = factor(flooded))) +
  geom_point(size = 3, position = position_dodge(0.2)) +
  geom_errorbar(aes(ymin = v_green_pct_lwr, ymax = v_green_pct_upr),
                width = 0.1, position = position_dodge(0.2)) +
  geom_line(aes(group = factor(flooded)), position = position_dodge(0.2)) +
  facet_wrap(~land_name, scales = "free_y") +
  scale_color_stata(name = "Affected?", labels = c("Yes", "No"), breaks = c(1,0)) +
  theme_minimal(base_size = 12) +
  labs(x = NULL, y = "Green Vote Share", title = "Primary Measure") +
  theme(legend.position = "top")
ggsave("~/Desktop/appendix_figure10.pdf", width = 7, height = 7)

df.fig11 <- 
  df.main %>% 
  group_by(date, flooded2, land_name) %>% 
  mutate(v2_green_pct_m = lm_robust(v_green_pct~1)$coefficients,
         v2_green_pct_lwr = lm_robust(v_green_pct~1)$conf.low,
         v2_green_pct_upr = lm_robust(v_green_pct~1)$conf.high) %>% 
  ungroup()
ggplot(df.fig11, aes(factor(date), v2_green_pct_m, color = factor(flooded2))) +
  geom_point(size = 3, position = position_dodge(0.2)) +
  geom_errorbar(aes(ymin = v2_green_pct_lwr, ymax = v2_green_pct_upr),
                width = 0.1, position = position_dodge(0.2)) +
  geom_line(aes(group = factor(flooded2)), position = position_dodge(0.2)) +
  facet_wrap(~land_name, scales = "free_y") +
  scale_color_stata(name = "Affected?", labels = c("Yes", "No"), breaks = c(1,0)) +
  theme_minimal(base_size = 12) +
  labs(x = NULL, y = "Green Vote Share", title = "Secondary Measure") +
  theme(legend.position = "top")
ggsave("~/Desktop/appendix_figure11.pdf", width = 7, height = 7)

# table 10 & figure 12  ---------------------------------------------------

## models ----
IDs.state <- unique(df.main$land_name)
tab10 <- vector("list", length = length(IDs.state)) %>% 
  `names<-`(IDs.state)
for (i in seq_along(tab10)) {
  if (names(tab10)[i]=="Sachsen") { 
    tab10[[i]][[1]] <- feols(
      data = filter(df.main, land_name == IDs.state[i]),
      fml = v_green_pct~factor(flooded)*factor(date)+pop_per_sqkm+net_out_per_thousand+old_share+income_mean+unemployed_rate+land_set_pct+land_agri_pct+factor(incumbent)*factor(date)+distance*factor(date), 
      fixef = "gen_id",
      cluster = "gen_id")
    tab10[[i]][[2]] <- NULL 
  }
  else {
    tab10[[i]][[1]] <- feols(
      data = filter(df.main, land_name == IDs.state[i]),
      fml = v_green_pct~factor(flooded)*factor(date)+pop_per_sqkm+net_out_per_thousand+old_share+income_mean+unemployed_rate+land_set_pct+land_agri_pct+distance*factor(date), 
      fixef = "gen_id",
      cluster = "gen_id")
    tab10[[i]][[2]] <- feols(
      data = filter(df.main, land_name == IDs.state[i]),
      fml = v_green_pct~factor(flooded2)*factor(date)+pop_per_sqkm+net_out_per_thousand+old_share+income_mean+unemployed_rate+land_set_pct+land_agri_pct+distance*factor(date), 
      fixef = "gen_id",
      cluster = "gen_id")
  }
}  

## table 10 ----
screenreg(c(map(tab10, 1), map(tab10, 2)[-4]), digits = 3, custom.header = list("Primary Measure"=1:4, "Secondary Measure"=5:7))

## figure 12 ----
df.fig12 <- 
  data.frame(
    state = rep(IDs.state, 2),
    term = rep(c("Primary", "Secondary"), each = 4),
    mu = c(0.003117315,-0.0004095530,0.0012558488,-0.005282131,0.006168334,0.0065201233,-0.0057605281,NA),
    se = c(0.0019963423,0.002237524,0.000951935,0.001427947,0.0036141077,0.002862208,0.001360638,NA)
  ) %>% 
  mutate(lwr=mu-1.96*se, upr=mu+1.96*se) %>% 
  mutate(state = factor(state, levels = sort(IDs.state)),
         term = factor(term, levels = IDs.measures))
ggplot(df.fig12) +
  geom_pointrange(aes(x="bhat", y=mu, ymin=lwr, ymax=upr, color=state, shape=state), position = position_dodge(0.7)) +
  facet_grid(rows = vars(term)) +
  geom_hline(yintercept = 0, color = "darkred", lty = "dotted") +
  coord_flip() +
  ggthemes::scale_color_stata() +
  labs(x = NULL, y = NULL, color = "State", shape = "State") +
  scale_y_continuous(breaks = scales::pretty_breaks(n=5)) +
  theme_minimal(base_size = 14) +
  theme(axis.ticks.y = element_blank(), axis.text.y = element_blank())
ggsave("~/Desktop/appendix_figure12.pdf", width = 9, height = 5)

# table 11 & figure 13 ----------------------------------------------------

## data ----
df.main %<>% mutate(nsevere = as.integer(flooded==1&severe==0),
                    nsevere2 = as.integer(flooded2==1&severe2==0))

## models ----
tab11 <- vector("list", length = length(IDs.state)) %>% 
  `names<-`(IDs.state)
for (i in seq_along(tab11)) {
  if (names(tab11)[i]=="Sachsen") { 
    tab11[[i]][[1]] <- tab11[[i]][[2]] <- NA 
  }
  if (names(tab11)[i]=="Bayern") { 
    tab11[[i]][[1]] <- feols( 
      data = filter(df.main, land_name == IDs.state[i]),
      fml = v_green_pct~factor(nsevere)*factor(date)+factor(severe)*factor(date)+pop_per_sqkm+net_out_per_thousand+old_share+income_mean+unemployed_rate+land_set_pct+land_agri_pct+distance*factor(date), 
      fixef = "gen_id",
      cluster = "gen_id") 
    tab11[[i]][[2]] <- NA 
  }
  else {
    tab11[[i]][[1]] <- feols(
      data = filter(df.main, land_name == IDs.state[i]),
      fml = v_green_pct~factor(nsevere)*factor(date)+factor(severe)*factor(date)+pop_per_sqkm+net_out_per_thousand+old_share+income_mean+unemployed_rate+land_set_pct+land_agri_pct+distance*factor(date), 
      fixef = "gen_id",
      cluster = "gen_id"
    )
    tab11[[i]][[2]] <- feols(
      data = filter(df.main, land_name == IDs.state[i]),
      fml = v_green_pct~factor(nsevere2)*factor(date)+factor(severe2)*factor(date)+pop_per_sqkm+net_out_per_thousand+old_share+income_mean+unemployed_rate+land_set_pct+land_agri_pct+distance*factor(date), 
      fixef = "gen_id",
      cluster = "gen_id")
  }
}

## table 11 ----
screenreg(c(map(tab11[1:3],1), map(tab11[1:2],2)), digits = 3, custom.header = list("Primary Measure"=1:3, "Secondary Measure"=4:5))
  
## figure 13 ----
df.fig13 <- 
  data.frame(
    state = rep(rep(IDs.state, each = 2), 2),
    term = rep(IDs.measures, each = 8),
    term1 = rep(c("Low","High"), 8),
    mu = c(0.006163193,0.000305194,-0.002017811,0.004977198,0.0021099842,-0.0030692848,NA,NA,0.004719248,0.001532154,0.0092261743,0.0040286889,NA,NA,NA,NA),
    se = c(0.0022890450,0.0026579141,0.002266321,0.003602483,0.001028285,0.001557863,NA,NA,0.0164006617,0.0167519667,0.004058736,0.003447607,NA,NA,NA,NA)
  ) %>% 
  mutate(lwr=mu-1.96*se, upr=mu+1.96*se) %>% 
  mutate(state = factor(state, levels = sort(IDs.state)),
         term = factor(term, levels = IDs.measures),
         term1 = factor(term1, levels = c("Low","High")))
ggplot(df.fig13) +
  geom_pointrange(aes(x="bhat", y=mu, ymin=lwr, ymax=upr, color=state, shape=term1), position = position_dodge(0.7)) +
  facet_grid(rows = vars(term)) +
  geom_hline(yintercept = 0, color = "darkred", lty = "dotted") +
  coord_flip() +
  ggthemes::scale_color_stata() +
  labs(x = NULL, y = NULL, color = "State", shape = "Severity") +
  scale_y_continuous(breaks = scales::pretty_breaks(n=5)) +
  theme_minimal(base_size = 14) +
  theme(axis.ticks.y = element_blank(), axis.text.y = element_blank())
ggsave("~/Desktop/appendix_figure13.pdf", width = 9, height = 5)

# figure 14 ---------------------------------------------------------------

df.s2.d <- 
  df.s2 %>% 
  filter(flooded2==1) %>% 
  group_by(gen_id) %>% 
  arrange(date) %>% 
  summarise(dv_green_pct = v_green_pct - lag(v_green_pct),
            mpop_per_sqkm = mean(pop_per_sqkm, na.rm=TRUE),
            severe2 = unique(severe2)) %>% 
  ungroup() %>% 
  filter(!is.na(dv_green_pct))
df.s2.d.t7 <- 
  df.s2.d %>% 
  arrange(desc(mpop_per_sqkm)) %>% 
  slice_head(n=7) %$% 
  filter(df.s2, gen_id %in% (.$gen_id)) %>% 
  select(gen_id, gen_name, date) %>% 
  filter(date==max(date)) %>% 
  distinct() %>% 
  left_join(df.s2.d %>% select(gen_id, dv_green_pct, mpop_per_sqkm))

ggplot(df.s2.d) +
  geom_point(aes(mpop_per_sqkm, dv_green_pct, color = factor(severe2), size = mpop_per_sqkm), 
             alpha = 0.7) +
  scale_color_manual(values = c("#006000","#ff4500"), 
                     breaks = c(0,1), 
                     labels = c("Low","High"), name = "Severity") +
  geom_text(data = df.s2.d.t7, 
            aes(mpop_per_sqkm, dv_green_pct, label = gen_name),
            nudge_y = 0.01) +
  labs(x = "Avg. Population Density", y = "Green Vote Change") +
  scale_x_continuous(breaks = scales::pretty_breaks(n=5)) +
  scale_y_continuous(labels = scales::percent) +
  guides(size = "none")
ggsave("~/Desktop/appendix_figure14.pdf", width = 7, height = 5)

# figure 15 ---------------------------------------------------------------

## data ----
df.eb <- read.csv("data_eb.csv", row.names = 1)
df.eb.m <- df.eb %>% 
  group_by(year, region_id1) %>% 
  summarise(
    climate_m = lm_robust(climate~1, weights = wts)$coefficients,
    climate_lwr = lm_robust(climate~1, weights = wts)$conf.low,
    climate_upr = lm_robust(climate~1, weights = wts)$conf.high)

## plot ----
ggplot(df.eb.m, aes(x = year, y = climate_m, color = region_id1, fill = region_id1)) +
  geom_point(data = df.eb, aes(y = climate), alpha = 0.2, size = 0.3, position = position_jitterdodge(dodge.width = 0.7, jitter.width = 0.7, jitter.height = 0.7)) +
  geom_pointrange(aes(ymin = climate_lwr, ymax = climate_upr), size = 0.6, position = position_dodge(width = 0.7)) +
  geom_errorbar(aes(ymin = climate_lwr, ymax = climate_upr), width = 0.3, position = position_dodge(width = 0.7)) +
  geom_line(aes(group = factor(region_id1)), size = 0.7, position = position_dodge(width = 0.7)) +
  scale_x_continuous(breaks = unique(df.eb.m$year)) +
  scale_y_continuous(breaks = 1:10, limits = c(1,10)) +
  scale_color_viridis_d(option = "H", begin = 0.1, end = 0.9) +
  scale_fill_viridis_d(option = "H", begin = 0.1, end = 0.9) +
  labs(title = "Public opinion on climate change in Germany",
       subtitle = "Eurobarometer survey, all states, 2009-2021",
       x = "", y = "Climate Concerns (1-10)",
       color = "Flooded?", fill = "Flooded?") +
  theme_minimal() + labs_pubr() +
  theme(plot.title = element_text(size = 16),
        plot.subtitle = element_text(size = 14),
        legend.position = "top")
ggsave("~/Desktop/appendix_figure15.pdf", width = 8, height = 6)

# table 12 & figure 16 ----------------------------------------------------

## data ----
df.h2 <- 
  df.pcpt0 %>% 
  group_by(land_id) %>% 
  mutate(flooded_land = max(flooded, na.rm = TRUE)) %>% 
  ungroup() %>% 
  filter(flooded_land == 1) 

## models ----
xx.h2 <- "+pop_per_sqkm+net_out_per_thousand+old_share+income_mean+unemployed_rate+land_set_pct+land_agri_pct+factor(incumbent)+distance*factor(year)"
fm.h2.a <- "climate_position~factor(flooded)*factor(year)"
fm.h2.b <- "climate_salience~factor(flooded)*factor(year)"
fm.h2.x.a <- paste0(fm.h2.a, xx.h2)
fm.h2.x.b <- paste0(fm.h2.b, xx.h2)
fm2.h2.a <- gsub("flooded", "flooded2", fm.h2.a)
fm2.h2.b <- gsub("flooded", "flooded2", fm.h2.b)
fm2.h2.x.a <- paste0(fm2.h2.a, xx.h2)
fm2.h2.x.b <- paste0(fm2.h2.b, xx.h2)

## table 12 ----
tab12.a <- map(
  list(fm.h2.a, fm.h2.x.a, fm2.h2.a, fm2.h2.x.a),
  ~fit_twfe(df = df.h2, fm = .x, fe = "wk_id", cl = "wk_id"))
tab12.b <- map(
  list(fm.h2.b, fm.h2.x.b, fm2.h2.b, fm2.h2.x.b),
  ~fit_twfe(df = df.h2, fm = .x, fe = "wk_id", cl = "wk_id"))
screenreg(c(tab12.a, tab12.b), digits = 3, custom.header = list("Climate Position"=1:4, "Issue Salience"=5:8))

## figure 16 ----
df.fig16 <- 
  list(tab12.a,tab12.b) %>% 
  `names<-`(c("h2.a","h2.b")) %>% 
  map_depth(2, ~{coeftable(.x) %>% data.frame}) %>% 
  map_depth(2, ~{rownames_to_column(.x, var = "term") %>% select(term, mean = Estimate, se = 3)}) %>% 
  map_depth(2, ~{filter(.x, term %in% c("factor(flooded)1:factor(year)2021", "factor(flooded2)1:factor(year)2021")) %>% 
      mutate(term = case_when(grepl("flooded2", term) ~ IDs.measures[2], TRUE~IDs.measures[1]))}) %>% 
  map(bind_rows) %>% 
  map(~add_column(.x, model = rep(IDs.models,2))) %>% 
  bind_rows(.id = "dv") %>% 
  mutate(lwr = mean-1.96*se, upr = mean+1.96*se) %>% 
  mutate(dv = ifelse(grepl("a", dv), "Position", "Salience"))
ggplot(df.fig16) +
  geom_pointrange(aes(x="bhat", y=mean, ymin=lwr, ymax=upr, color=model, shape=model), position = position_dodge(0.7)) +
  facet_grid(rows = vars(term), cols = vars(dv), scales = "free_x") +
  geom_hline(yintercept = 0, color = "darkred", lty = "dotted") +
  coord_flip() +
  scale_color_stata() +
  labs(x = NULL, y = NULL, color = "Model", shape = "Model", title = "Climate Concern Effect") +
  scale_y_continuous(breaks = scales::pretty_breaks(n=5)) +
  theme(axis.ticks.y = element_blank(), axis.text.y = element_blank())
ggsave("~/Desktop/appendix_figure16.pdf", width = 7, height = 4) 

# table 13 & figure 17 ----------------------------------------------------

## data ----
df.h3a <- readRDS("data_turnout_main.RDS") %>% 
  filter(gen_id %in% df.main$gen_id) %>% 
  left_join(df.main %>% filter(date == max(date)) %>% 
              select(gen_id, starts_with("flood"), v_green_pct))%>% 
  left_join(df.main0 %>% filter(date == max(date)) %>% 
              select(gen_id, starts_with("flood"), v_green_pct))

## table 13 ----
tab13 <- 
  rbind.data.frame(
    df.h3a %>% nest(-flooded) %>% 
      mutate(corr = map(data, ~cor.test(.$turnout, .x$v_green_pct)) %>% map(., broom::tidy)) %>% 
      unnest(corr) %>% select(-data) %>% add_column(measure = "primary"),
    df.h3a %>% nest(-flooded2) %>% 
      mutate(corr = map(data, ~cor.test(.$turnout, .x$v_green_pct)) %>% map(., broom::tidy)) %>% 
      unnest(corr) %>% select(-data) %>% rename(flooded = flooded2) %>% add_column(measure = "secondary")
  )
tab13[,c("flooded","estimate","statistic","p.value","measure")] %>% knitr::kable()

## figure 17 ----
tab13 %>% 
  mutate(flooded = factor(flooded, levels = c(1,0), labels = c("Yes","No")),
         measure = str_to_title(measure)) %>% 
  ggplot() +
  geom_col(aes(x = flooded, y = estimate, color = flooded, fill = flooded)) +
  geom_errorbar(aes(x = flooded, y = estimate, ymin = conf.low, ymax = conf.high), width = 0.1) +
  facet_grid(rows = vars(measure)) +
  geom_hline(yintercept = 0, color = "darkred", lty = "dotted") +
  coord_flip() + 
  scale_color_manual(values = c("No" = "#1a476f", "Yes" = "#90353b")) +
  scale_fill_manual(values = c("No" = "#1a476f", "Yes" = "#90353b")) +
  labs(x = NULL, y = NULL,  color = "Flooded", fill = "Flooded", title = "Mobilisation Effect") +
  scale_y_continuous(breaks = scales::pretty_breaks(n=5)) +
  theme(axis.ticks.y = element_blank(), axis.text.y = element_blank())
ggsave("~/Desktop/appendix_figure17.pdf", width = 7, height = 4)  

# table 14 ----------------------------------------------------------------

## data ----
df.med.turnout <- rbind.data.frame(
  readRDS("data_turnout_main.RDS") %>% 
    select(gen_id, date, turnout, date),
  readRDS("data_turnout_ts.RDS") %>% 
    filter(date==max(date)) %>% 
    mutate(turnout = turnout/100) %>% 
    select(gen_id, date, turnout)
)
df.med.turnout.fd <- 
  df.med.turnout %>% 
  mutate(date=as.Date(date), gen_id = as.integer(gen_id)) %>% 
  group_by(gen_id) %>% 
  summarise(turnout = turnout - lag(turnout)) %>% 
  ungroup() %>% 
  na.omit()
df.med.fd <- 
  df.main0 %>% 
  mutate(date=as.Date(date), gen_id = as.integer(gen_id)) %>% 
  group_by(gen_id) %>% 
  summarise(across(.cols = c(v_green_pct, pop_per_sqkm:incumbent),
                   .fns = ~{.x - lag(.x)})) %>% 
  ungroup() %>% 
  na.omit()
df.med.turnouts <- 
  left_join(df.med.fd, df.med.turnout.fd, by = "gen_id") %>% 
  left_join(
    df.main0 %>% dplyr::select(gen_id, flooded, flooded2) %>% distinct(),
    by = "gen_id"
  ) %>% 
  mutate(
    flooded=factor(flooded), 
    flooded2=factor(flooded2))
df.med.turnouts.s4 <- 
  df.med.turnouts %>% 
  filter(gen_id %in% df.main$gen_id)

## models ----
pacman::p_load(mediation) # avoid ^namespace conflict (MASS::select vs dplyr::select)
set.seed(1)

tab14.med <- lm(turnout~flooded+pop_per_sqkm+net_out_per_thousand+old_share+income_mean+unemployed_rate+land_set_pct+land_agri_pct+factor(incumbent), data = df.med.turnouts.s4)
tab14.out <- lm(v_green_pct~turnout+flooded+pop_per_sqkm+net_out_per_thousand+old_share+income_mean+unemployed_rate+land_set_pct+land_agri_pct+factor(incumbent), data = df.med.turnouts.s4)
tab14.1 <- mediate(tab14.med, tab14.out, treat = "flooded", mediator = "turnout", robustSE = TRUE)
tab14.med2 <- lm(turnout~flooded2+pop_per_sqkm+net_out_per_thousand+old_share+income_mean+unemployed_rate+land_set_pct+land_agri_pct+factor(incumbent), data = df.med.turnouts.s4)
tab14.out2 <- lm(v_green_pct~turnout+flooded2+pop_per_sqkm+net_out_per_thousand+old_share+income_mean+unemployed_rate+land_set_pct+land_agri_pct+factor(incumbent), data = df.med.turnouts.s4)
tab14.2 <- mediate(tab14.med2, tab14.out2, treat = "flooded2", mediator = "turnout", robustSE = TRUE)
tab14 <- list(tab14.1,tab14.2)

pacman::p_unload(mediation)
select <- dplyr::select # enforce namespace clarity

## summaries ----
summary(tab14.1)
summary(tab14.2)

# table 15 & figure 18 ----------------------------------------------------

## data ----
df.h3b <- read.csv("data_gles_individual.csv", row.names = 1) %>% 
  mutate(v_change = v_2021-v_2017) %>% 
  group_by(land_id) %>% 
  mutate(flooded_land = max(flooded, na.rm = TRUE)) %>% 
  ungroup() %>% 
  filter(flooded_land == 1)

## models ----
xx.h3b <- "+factor(gender)+age+factor(married)+factor(immigrant)+educ+factor(empl)+hhinc+urban+relig+lr"
fm.h3b <- "v_change~factor(flooded)"
fm.h3b.x <- paste0(fm.h3b, xx.h3b)
fm2.h3b <- "v_change~factor(flooded2)"
fm2.h3b.x <- paste0(fm2.h3b, xx.h3b)

## table 15 ----
tab15 <- map(
  list(fm.h3b, fm.h3b.x, fm2.h3b, fm2.h3b.x),
  ~fit_twfe2(df = df.h3b, fml = .x, cl = "wk_id", wt = df.h3b$w_ow))
screenreg(tab15, digits = 3)

## figure 18 ----
df.fig18 <- 
  tab15 %>% 
  map(~{coeftable(.x) %>% as.data.frame}) %>% 
  map(~{rownames_to_column(.x, var = "term") %>% select(term, mean = Estimate, se = `Std. Error`)}) %>% 
  map(~{filter(.x, term %in% c("factor(flooded)1", "factor(flooded2)1")) %>% 
      mutate(term = case_when(grepl("flooded2", term) ~ IDs.measures[2], TRUE~IDs.measures[1]))}) %>% 
  `names<-`(rep(IDs.models,2)) %>% 
  bind_rows(.id = "model") %>% 
  mutate(lwr = mean-1.96*se, upr = mean+1.96*se)
ggplot(df.fig18) +
  geom_pointrange(aes(x="bhat", y=mean, ymin=lwr, ymax=upr, color=model, shape = model), position = position_dodge(0.7)) +
  facet_grid(rows = vars(term)) +
  geom_hline(yintercept = 0, color = "darkred", lty = "dotted") +
  coord_flip() +
  scale_color_stata() +
  labs(x = NULL, y = NULL, color = "Model", shape = "Model", title = "Persuasion Effect") +
  scale_y_continuous(limits = c(-0.1,0.1), breaks = scales::pretty_breaks(n=5)) +
  theme(axis.ticks.y = element_blank(), axis.text.y = element_blank())
ggsave("~/Desktop/appendix_figure18.pdf", width = 7, height = 4)  

# tables 16 & 17 ----------------------------------------------------------

## data ----
df.h3b %<>%   
  filter(land_id %in% c(5, 7, 9, 14)) %>% 
  mutate(v1_2017 = relevel(factor(v1_2017), ref = "Others"))

## models ----
fm.swing <- "v_2021~factor(v1_2017)*factor(flooded)+factor(gender)+age+factor(married)+factor(immigrant)+educ+factor(empl)+hhinc+urban+relig+lr"
fm2.swing <- gsub("flooded", "flooded2", fm.swing)

## table 16 ----
tab16 <- map(
  list(fm.swing, fm2.swing),
  ~fit_twfe2(df = df.h3b, fml = .x, fe = "wk_id", cl = "wk_id", wt = df.h3b$w_ow)) 
screenreg(tab16, digits = 3)

## table 17 ----
tab17 <- map(
  list(fm.swing, fm2.swing),
  ~fit_twfe2_logit(df = df.h3b, fml = .x, fe = "wk_id", cl = "wk_id", wt = df.h3b$w_ow)) 
screenreg(tab17, digits = 3)

# cleanup -----------------------------------------------------------------

theme_set(theme.default)
pacman::p_unload("all")
rm(list = ls())
